Metrics for assessing community assembly processes

source('../env.R')

Skipping install of 'clootl' from a github remote, the SHA1 (2ed1650b) has not changed since last install.
  Use `force = TRUE` to force installation

community_data = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'community_assembly_metrics_using_relative_abundance.csv'))

Rows: 308 Columns: 10── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
dbl (10): mntd_standard, mntd_actual, mass_fdiv_standard, mass_fdiv_actual, beak_width_fdiv_standard, beak_width_fdiv_actual, hwi_fdiv_standard, hwi_fdiv_act...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(community_data)

colnames(community_data)

 [1] "mntd_standard"            "mntd_actual"              "mass_fdiv_standard"       "mass_fdiv_actual"         "beak_width_fdiv_standard"
 [6] "beak_width_fdiv_actual"   "hwi_fdiv_standard"        "hwi_fdiv_actual"          "city_id"                  "urban_pool_size"         

min(community_data$mntd_standard)

[1] -2.33692

max(community_data$mntd_standard)

[1] 2.328448

min(community_data$beak_width_fdiv_standard)

[1] -2.685152

max(community_data$beak_width_fdiv_standard)

[1] 1.931681

min(community_data$hwi_fdiv_standard)

[1] -2.200336

max(community_data$hwi_fdiv_standard)

[1] 2.333383

min(community_data$mass_fdiv_standard)

[1] -2.377212

max(community_data$mass_fdiv_standard)

[1] 2.1073

Join on realms

city_to_realm = read_csv(filename(CITY_DATA_OUTPUT_DIR, 'realms.csv'))

Rows: 337 Columns: 2── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): core_realm
dbl (1): city_id
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

community_data_with_realm = left_join(community_data, city_to_realm)

Joining with `by = join_by(city_id)`

Cities as points

city_points = st_centroid(read_sf(filename(CITY_DATA_OUTPUT_DIR, 'city_selection.shp'))) %>% left_join(community_data_with_realm)

Warning: st_centroid assumes attributes are constant over geometriesWarning: st_centroid does not give correct centroids for longitude/latitude dataJoining with `by = join_by(city_id)`

city_points_coords = st_coordinates(city_points)
city_points$latitude = city_points_coords[,1]
city_points$longitude = city_points_coords[,2]

world_map = read_country_boundaries()

Load community data, and create long format version

communities = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'communities_for_analysis.csv'))

Rows: 2428 Columns: 7── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (5): city_name, ebird_species_name, seasonal, presence, origin
dbl (2): city_id, relative_abundance_proxy
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

communities

community_summary = communities %>% group_by(city_id) %>% summarise(regional_pool_size = n(), urban_pool_size = sum(relative_abundance_proxy > 0))
community_summary

Load trait data

traits = read_csv(filename(TAXONOMY_OUTPUT_DIR, 'traits_ebird.csv'))

Rows: 332 Columns: 10── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (5): ebird_species_name, habitat, trophic_level, trophic_niche, primary_lifestyle
dbl (5): beak_width, hwi, mass, habitat_density, migration
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(traits)

Load realm geo

resolve = read_resolve()
head(resolve)

Simple feature collection with 6 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -162.1547 ymin: -69.55876 xmax: 158.6167 ymax: 61.53428
Geodetic CRS:  WGS 84

Load spatial var

spatial_var = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'spatial_var.csv')) %>% filter(city_id %in% community_summary$city_id)

Rows: 337 Columns: 3── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
dbl (3): city_id, NMDS1, NMDS2
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

spatial_var

Summary metrics by Realm

test_required_values = function(name, df) {
  cat(paste(
    test_value_wilcox(paste(name, 'MNTD'), df$mntd_standard),
    test_value_wilcox(paste(name, 'Beak Width FDiv'), df$beak_width_fdiv_standard),
    test_value_wilcox(paste(name, 'HWI FDiv'), df$hwi_fdiv_standard),
    test_value_wilcox(paste(name, 'Mass FDiv'), df$mass_fdiv_standard),
    paste('N', nrow(df)),
    sep = "\n"))
}

test_required_values('Global', community_data_with_realm)

Global MNTD median -0.36 ***
Global Beak Width FDiv median 0.02 
Global HWI FDiv median 0.39 **
Global Mass FDiv median 0.29 ***
N 308

unique(community_data_with_realm$core_realm)

[1] "Nearctic"    "Neotropic"   "Palearctic"  "Afrotropic"  "Indomalayan" "Australasia"

test_required_values('Nearctic', community_data_with_realm[community_data_with_realm$core_realm == 'Nearctic',])

Nearctic MNTD median 0.67 *
Nearctic Beak Width FDiv median 0.29 
Nearctic HWI FDiv median -0.8 ***
Nearctic Mass FDiv median -0.26 
N 46

test_required_values('Neotropic', community_data_with_realm[community_data_with_realm$core_realm == 'Neotropic',])

Neotropic MNTD median 0.03 
Neotropic Beak Width FDiv median -0.44 ***
Neotropic HWI FDiv median -0.31 
Neotropic Mass FDiv median 0.33 *
N 64

test_required_values('Palearctic', community_data_with_realm[community_data_with_realm$core_realm == 'Palearctic',])

Palearctic MNTD median 0.13 
Palearctic Beak Width FDiv median 1.25 ***
Palearctic HWI FDiv median -0.39 
Palearctic Mass FDiv median 0.01 
N 72

test_required_values('Afrotropic', community_data_with_realm[community_data_with_realm$core_realm == 'Afrotropic',])

Afrotropic MNTD median -1.28 *
Afrotropic Beak Width FDiv median -0.56 
Afrotropic HWI FDiv median 0.15 
Afrotropic Mass FDiv median -0.95 
N 9

test_required_values('Indomalayan', community_data_with_realm[community_data_with_realm$core_realm == 'Indomalayan',])

Indomalayan MNTD median -0.64 ***
Indomalayan Beak Width FDiv median -0.68 ***
Indomalayan HWI FDiv median 1.11 ***
Indomalayan Mass FDiv median 0.83 ***
N 111

test_required_values('Australasia', community_data_with_realm[community_data_with_realm$core_realm == 'Australasia',])

Australasia MNTD median -1.39 
Australasia Beak Width FDiv median -0.75 
Australasia HWI FDiv median 0.77 
Australasia Mass FDiv median -0.96 
N 6

Summary metrics by invasive species

cities_with_introduced_species = communities %>% filter(origin == 'Introduced') %>% select(city_id) %>% distinct()

cities_with_no_introduced_species = communities %>% filter(!(city_id %in% cities_with_introduced_species$city_id)) %>% select(city_id) %>% distinct()

cities_with_introduced_species$introduced_species = TRUE
cities_with_no_introduced_species$introduced_species = FALSE

community_data_with_realm_with_introduced = community_data_with_realm %>% left_join(rbind(cities_with_introduced_species, cities_with_no_introduced_species))

Joining with `by = join_by(city_id)`

community_data_with_realm_with_introduced

test_required_values('With Introduced', community_data_with_realm_with_introduced[community_data_with_realm_with_introduced$introduced_species,])

With Introduced MNTD median -0.03 
With Introduced Beak Width FDiv median 0.11 
With Introduced HWI FDiv median -0.36 
With Introduced Mass FDiv median 0.01 
N 189

test_required_values('Without Introduced', community_data_with_realm_with_introduced[!community_data_with_realm_with_introduced$introduced_species,])

Without Introduced MNTD median -0.53 ***
Without Introduced Beak Width FDiv median -0.28 *
Without Introduced HWI FDiv median 1.04 ***
Without Introduced Mass FDiv median 0.72 ***
N 119

What families exist in which realms?

communities %>% 
  left_join(city_to_realm) %>% 
  mutate(family = gsub( " .*$", "", ebird_species_name)) %>%
  dplyr::select(family, core_realm) %>%
  distinct() %>%
  arrange(core_realm)

Joining with `by = join_by(city_id)`

Total urban families

communities %>% 
  mutate(family = gsub( " .*$", "", ebird_species_name)) %>%
  dplyr::select(family) %>%
  distinct() %>%
  arrange()

of which urban

communities %>% 
  filter(relative_abundance_proxy > 0) %>%
  mutate(family = gsub( " .*$", "", ebird_species_name)) %>%
  dplyr::select(family) %>%
  distinct() %>%
  arrange()

Summary metrics by introduced species

communities = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'communities_for_analysis.csv'))

Rows: 2428 Columns: 7── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (5): city_name, ebird_species_name, seasonal, presence, origin
dbl (2): city_id, relative_abundance_proxy
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

city_introduced_species = communities %>% group_by(city_id) %>% summarise(number_of_species = n()) %>% left_join(
  communities %>% group_by(city_id) %>% filter(origin == 'Introduced') %>% summarise(number_of_introduced_species = n())
) %>% replace_na(list(number_of_introduced_species = 0))

Joining with `by = join_by(city_id)`

community_data_with_introductions = left_join(community_data, city_introduced_species)

Joining with `by = join_by(city_id)`

community_data_with_introductions$has_introduced_species = community_data_with_introductions$number_of_introduced_species > 0
community_data_with_introductions

communities %>% 
  filter(origin == 'Introduced') %>%
  dplyr::select(ebird_species_name) %>%
  group_by(ebird_species_name) %>%
  summarise(total_cities = n()) %>%
  arrange(desc(total_cities))

community_data_with_introductions[,c('mntd_standard', 'has_introduced_species')]

community_data_with_introductions %>% group_by(has_introduced_species) %>% summarise(
  total_cities = n(), 
  
  mean_mntd_std = mean(mntd_standard, na.rm = T),
  median_mntd_std = median(mntd_standard, na.rm = T),
  sd_mntd_std = sd(mntd_standard, na.rm = T),
  
  mean_mass_fdiv_std = mean(mass_fdiv_standard, na.rm = T),
  median_mass_fdiv_std = median(mass_fdiv_standard, na.rm = T),
  sd_mass_fdiv_std = sd(mass_fdiv_standard, na.rm = T),
  
  mean_gape_width_fdiv_std = mean(beak_width_fdiv_standard, na.rm = T),
  median_gape_width_fdiv_std = median(beak_width_fdiv_standard, na.rm = T),
  sd_gape_width_fdiv_std = sd(beak_width_fdiv_standard, na.rm = T),
  
  mean_handwing_index_fdiv_std = mean(hwi_fdiv_standard, na.rm = T),
  median_handwing_index_fdiv_std = median(hwi_fdiv_standard, na.rm = T),
  sd_handwing_index_fdiv_std = sd(hwi_fdiv_standard, na.rm = T)
)

MNTD

ggplot(community_data_with_introductions, aes(x = has_introduced_species, y = mntd_standard)) + geom_boxplot()

wilcox.test(mntd_standard ~ has_introduced_species, community_data_with_introductions, na.action = 'na.omit')

    Wilcoxon rank sum test with continuity correction

data:  mntd_standard by has_introduced_species
W = 7925, p-value = 0.00001285
alternative hypothesis: true location shift is not equal to 0

There is a significant difference between the response of cities with introduced species (0.53±0.27) and those without (0.47±0.19) (p-value = 0.02).

Mass FDiv

ggplot(community_data_with_introductions, aes(x = has_introduced_species, y = mass_fdiv_standard)) + geom_boxplot()

wilcox.test(mass_fdiv_standard ~ has_introduced_species, community_data_with_introductions, na.action = 'na.omit')

    Wilcoxon rank sum test with continuity correction

data:  mass_fdiv_standard by has_introduced_species
W = 15028, p-value = 0.0000006706
alternative hypothesis: true location shift is not equal to 0

There is a significant difference between the response of cities with introduced species (0.57±0.27) and those without (0.73±0.24) (p < 0.0001)

Beak Gape FDiv

ggplot(community_data_with_introductions, aes(x = has_introduced_species, y = beak_width_fdiv_standard)) + geom_boxplot()

wilcox.test(beak_width_fdiv_standard ~ has_introduced_species, community_data_with_introductions, na.action = 'na.omit')

    Wilcoxon rank sum test with continuity correction

data:  beak_width_fdiv_standard by has_introduced_species
W = 8662, p-value = 0.0006884
alternative hypothesis: true location shift is not equal to 0

There is NOT a significant difference between the response of cities with introduced species (0.61±0.30) and those without (0.56±0.27)

HWI FDiv

ggplot(community_data_with_introductions, aes(x = has_introduced_species, y = hwi_fdiv_standard)) + geom_boxplot()

wilcox.test(hwi_fdiv_standard ~ has_introduced_species, community_data_with_introductions, na.action = 'na.omit')

    Wilcoxon rank sum test with continuity correction

data:  hwi_fdiv_standard by has_introduced_species
W = 17606, p-value < 0.00000000000000022
alternative hypothesis: true location shift is not equal to 0

There is a significant difference between the response of cities with introduced species (0.49±0.30) and those without (0.79±0.21) (p < 0.0001)

What proportion of cities in each realm have introduced species

community_data_with_introductions %>% left_join(city_to_realm) %>%
  group_by(core_realm) %>%
  summarise(
    total_cities = n(), 
    total_cities_with_introduced = sum(has_introduced_species), 
    total_cities_without_introduced = n() - sum(has_introduced_species)) %>%
  arrange(core_realm)

Joining with `by = join_by(city_id)`

Are any introduced species not also present in a city?

communities %>% 
  filter(origin == 'Introduced') %>% 
  filter(relative_abundance_proxy < 0.1)

What’s the average relative abundance of introduced species compared to native

communities %>% 
  group_by(origin) %>% 
  summarise(average_relative_abundance = mean(relative_abundance_proxy))

communities %>% 
  group_by(origin) %>% 
  filter(relative_abundance_proxy > 0) %>%
  summarise(average_relative_abundance = mean(relative_abundance_proxy))

communities %>% 
  group_by(origin) %>% 
  summarise(average_relative_abundance = mean(relative_abundance_proxy))

Create analysis data frame

geography = read_csv(filename(CITY_DATA_OUTPUT_DIR, 'geography.csv'))

Rows: 342 Columns: 26── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
dbl (26): city_id, city_avg_ndvi, city_avg_elevation, city_avg_temp, city_avg_min_monthly_temp, city_avg_max_monthly_temp, city_avg_monthly_temp, city_avg_ra...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

names(geography)

 [1] "city_id"                       "city_avg_ndvi"                 "city_avg_elevation"            "city_avg_temp"                
 [5] "city_avg_min_monthly_temp"     "city_avg_max_monthly_temp"     "city_avg_monthly_temp"         "city_avg_rainfall"            
 [9] "city_avg_max_monthly_rainfall" "city_avg_min_monthly_rainfall" "city_avg_soil_moisture"        "city_max_elev"                
[13] "city_min_elev"                 "city_elev_range"               "region_20km_avg_ndvi"          "region_20km_avg_elevation"    
[17] "region_20km_avg_soil_moisture" "region_20km_max_elev"          "region_20km_min_elev"          "region_20km_elev_range"       
[21] "region_50km_avg_ndvi"          "region_50km_avg_elevation"     "region_50km_avg_soil_moisture" "region_50km_max_elev"         
[25] "region_50km_min_elev"          "region_50km_elev_range"       

analysis_data = community_data_with_realm[,c('city_id', 'mntd_standard', 'mass_fdiv_standard', 'beak_width_fdiv_standard', 'hwi_fdiv_standard', 'core_realm')] %>% 
  left_join(city_points[,c('city_id', 'latitude', 'longitude')]) %>%
  left_join(community_data_with_introductions[,c('city_id', 'has_introduced_species')]) %>%
  left_join(geography) %>%
  left_join(spatial_var)

Joining with `by = join_by(city_id)`Joining with `by = join_by(city_id)`Joining with `by = join_by(city_id)`Joining with `by = join_by(city_id)`

analysis_data$abs_latitude = abs(analysis_data$latitude)
analysis_data$core_realm = factor(analysis_data$core_realm, levels = c('Palearctic', 'Nearctic', 'Neotropic', 'Afrotropic', 'Indomalayan', 'Australasia', 'Oceania'))
analysis_data$has_introduced_species = factor(analysis_data$has_introduced_species, level = c('FALSE', 'TRUE'), labels = c('No introduced species', 'Introduced species'))

model_data = function(df, dependant_var) {
  df[,c(dependant_var, 'core_realm', 'abs_latitude', 'longitude', 'has_introduced_species', 'city_avg_ndvi', 'city_avg_elevation', 'city_avg_temp', 'city_avg_min_monthly_temp', 'city_avg_max_monthly_temp', 'city_avg_monthly_temp', 'city_avg_rainfall', 'city_avg_max_monthly_rainfall', 'city_avg_min_monthly_rainfall', 'city_avg_soil_moisture', 'city_max_elev', 'city_min_elev', 'city_elev_range', 'region_20km_avg_ndvi', 'region_20km_avg_elevation', 'region_20km_avg_soil_moisture', 'region_20km_max_elev', 'region_20km_min_elev', 'region_20km_elev_range', 'region_50km_avg_ndvi', 'region_50km_avg_elevation', 'region_50km_avg_soil_moisture', 'region_50km_max_elev', 'region_50km_min_elev', 'region_50km_elev_range')]
}
model_data(analysis_data, 'mntd_standard')

names(analysis_data)

 [1] "city_id"                       "mntd_standard"                 "mass_fdiv_standard"            "beak_width_fdiv_standard"     
 [5] "hwi_fdiv_standard"             "core_realm"                    "latitude"                      "longitude"                    
 [9] "geometry"                      "has_introduced_species"        "city_avg_ndvi"                 "city_avg_elevation"           
[13] "city_avg_temp"                 "city_avg_min_monthly_temp"     "city_avg_max_monthly_temp"     "city_avg_monthly_temp"        
[17] "city_avg_rainfall"             "city_avg_max_monthly_rainfall" "city_avg_min_monthly_rainfall" "city_avg_soil_moisture"       
[21] "city_max_elev"                 "city_min_elev"                 "city_elev_range"               "region_20km_avg_ndvi"         
[25] "region_20km_avg_elevation"     "region_20km_avg_soil_moisture" "region_20km_max_elev"          "region_20km_min_elev"         
[29] "region_20km_elev_range"        "region_50km_avg_ndvi"          "region_50km_avg_elevation"     "region_50km_avg_soil_moisture"
[33] "region_50km_max_elev"          "region_50km_min_elev"          "region_50km_elev_range"        "NMDS1"                        
[37] "NMDS2"                         "abs_latitude"                 

NMDS Spatial Helpers

analysis_data_nmds_coords = analysis_data[,c('NMDS1', 'NMDS2')]
coordinates(analysis_data_nmds_coords)  = ~ NMDS1 + NMDS2

analysis_data_nmds_nearneigh <- knearneigh(analysis_data_nmds_coords)

Warning: knearneigh: identical points foundWarning: knearneigh: kd_tree not available for identical points

analysis_data_nmds_neighbours <- knn2nb(analysis_data_nmds_nearneigh)

Polygons around realms in NMDS plot

cities_to_realms_nmds = read_csv(filename(CITY_DATA_OUTPUT_DIR, 'realms.csv')) %>% left_join(analysis_data) %>% filter(!is.na(NMDS1))

Rows: 337 Columns: 2── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): core_realm
dbl (1): city_id
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Joining with `by = join_by(city_id, core_realm)`

unique(cities_to_realms_nmds$core_realm)

[1] "Nearctic"    "Neotropic"   "Palearctic"  "Afrotropic"  "Indomalayan" "Australasia"

realm_nmds_neartic_polygon = cities_to_realms_nmds %>% filter(core_realm == 'Nearctic') %>% slice(chull(NMDS1, NMDS2))
realm_nmds_neotropic_polygon = cities_to_realms_nmds %>% filter(core_realm == 'Neotropic') %>% slice(chull(NMDS1, NMDS2))
realm_nmds_palearctic_polygon = cities_to_realms_nmds %>% filter(core_realm == 'Palearctic') %>% slice(chull(NMDS1, NMDS2))
realm_nmds_afrotropic_polygon = cities_to_realms_nmds %>% filter(core_realm == 'Afrotropic') %>% slice(chull(NMDS1, NMDS2))
realm_nmds_indomalayan_polygon = cities_to_realms_nmds %>% filter(core_realm == 'Indomalayan') %>% slice(chull(NMDS1, NMDS2))
realm_nmds_australasia_polygon = cities_to_realms_nmds %>% filter(core_realm == 'Australasia') %>% slice(chull(NMDS1, NMDS2))

polygon_line_type = 'dashed'
polygon_linewidth = 0.4

with_realms_nmds = function(g) {
  g + 
    geom_polygon(data = realm_nmds_neartic_polygon, linewidth = polygon_linewidth, linetype = polygon_line_type, alpha = 0) +
    geom_polygon(data = realm_nmds_neotropic_polygon, linewidth = polygon_linewidth, linetype = polygon_line_type, alpha = 0) +
    geom_polygon(data = realm_nmds_palearctic_polygon, linewidth = polygon_linewidth, linetype = polygon_line_type, alpha = 0) +
    geom_polygon(data = realm_nmds_afrotropic_polygon, linewidth = polygon_linewidth, linetype = polygon_line_type, alpha = 0) +
    geom_polygon(data = realm_nmds_indomalayan_polygon, linewidth = polygon_linewidth, linetype = polygon_line_type, alpha = 0) +
    geom_polygon(data = realm_nmds_australasia_polygon, linewidth = polygon_linewidth, linetype = polygon_line_type, alpha = 0)
}

Lat Long Spatial Helpers

analysis_data_latlong_coords = analysis_data[,c('longitude', 'latitude')]
coordinates(analysis_data_latlong_coords)  = ~ longitude + latitude

analysis_data_coords_nearneigh <- knearneigh(analysis_data_latlong_coords, longlat = TRUE)
analysis_data_neighbours <- knn2nb(analysis_data_coords_nearneigh)

Polygons around realms in lat long plot

cities_to_realms_latlong = read_csv(filename(CITY_DATA_OUTPUT_DIR, 'realms.csv')) %>% left_join(analysis_data) %>% filter(!is.na(latitude))

Rows: 337 Columns: 2── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): core_realm
dbl (1): city_id
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Joining with `by = join_by(city_id, core_realm)`

unique(cities_to_realms_latlong$core_realm)

[1] "Nearctic"    "Neotropic"   "Palearctic"  "Afrotropic"  "Indomalayan" "Australasia"

realm_latlong_neartic_polygon = cities_to_realms_latlong %>% filter(core_realm == 'Nearctic') %>% slice(chull(latitude, longitude))
realm_latlong_neotropic_polygon = cities_to_realms_latlong %>% filter(core_realm == 'Neotropic') %>% slice(chull(latitude, longitude))
realm_latlong_palearctic_polygon = cities_to_realms_latlong %>% filter(core_realm == 'Palearctic') %>% slice(chull(latitude, longitude))
realm_latlong_afrotropic_polygon = cities_to_realms_latlong %>% filter(core_realm == 'Afrotropic') %>% slice(chull(latitude, longitude))
realm_latlong_indomalayan_polygon = cities_to_realms_latlong %>% filter(core_realm == 'Indomalayan') %>% slice(chull(latitude, longitude))
realm_latlong_australasia_polygon = cities_to_realms_latlong %>% filter(core_realm == 'Australasia') %>% slice(chull(latitude, longitude))

with_realms_latlong = function(g) {
  g + 
    geom_polygon(data = realm_latlong_neartic_polygon, linewidth = polygon_linewidth, linetype = polygon_line_type, alpha = 0) +
    geom_polygon(data = realm_latlong_neotropic_polygon, linewidth = polygon_linewidth, linetype = polygon_line_type, alpha = 0) +
    geom_polygon(data = realm_latlong_palearctic_polygon, linewidth = polygon_linewidth, linetype = polygon_line_type, alpha = 0) +
    geom_polygon(data = realm_latlong_afrotropic_polygon, linewidth = polygon_linewidth, linetype = polygon_line_type, alpha = 0) +
    geom_polygon(data = realm_latlong_indomalayan_polygon, linewidth = polygon_linewidth, linetype = polygon_line_type, alpha = 0) +
    geom_polygon(data = realm_latlong_australasia_polygon, linewidth = polygon_linewidth, linetype = polygon_line_type, alpha = 0)
}

Check for spatial autocorrelation

MNTD

Lat/Long

with_realms_latlong(ggplot(analysis_data, aes(x = latitude, y = longitude, colour = mntd_standard)) + geom_point() + standardised_colours_scale + labs(colour = "Standardised response"))

moran.test(analysis_data$mntd_standard, nb2listw(analysis_data_neighbours))

    Moran I test under randomisation

data:  analysis_data$mntd_standard  
weights: nb2listw(analysis_data_neighbours)    

Moran I statistic standard deviate = 7.3361, p-value = 0.0000000000001099
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.529289710      -0.003257329       0.005269629 

NMDS

with_realms_nmds(ggplot(analysis_data, aes(x = NMDS1, y = NMDS2, colour = mntd_standard)) + geom_point() + standardised_colours_scale + labs(colour = "Standardised response"))

moran.test(analysis_data$mntd_standard, nb2listw(analysis_data_nmds_neighbours))

    Moran I test under randomisation

data:  analysis_data$mntd_standard  
weights: nb2listw(analysis_data_nmds_neighbours)    

Moran I statistic standard deviate = 6.5062, p-value = 0.00000000003854
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.444629830      -0.003257329       0.004738970 

FDiv Beak Width

Lat/Long

with_realms_latlong(ggplot(analysis_data, aes(x = latitude, y = longitude, colour = beak_width_fdiv_standard)) + geom_point() + standardised_colours_scale + labs(colour = "Standardised response"))

moran.test(analysis_data$beak_width_fdiv_standard, nb2listw(analysis_data_neighbours))

    Moran I test under randomisation

data:  analysis_data$beak_width_fdiv_standard  
weights: nb2listw(analysis_data_neighbours)    

Moran I statistic standard deviate = 7.754, p-value = 0.000000000000004451
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.560232185      -0.003257329       0.005281006 

NMDS

with_realms_nmds(ggplot(analysis_data, aes(x = NMDS1, y = NMDS2, colour = beak_width_fdiv_standard)) + geom_point() + standardised_colours_scale + labs(colour = "Standardised response"))

moran.test(analysis_data$beak_width_fdiv_standard, nb2listw(analysis_data_nmds_neighbours))

    Moran I test under randomisation

data:  analysis_data$beak_width_fdiv_standard  
weights: nb2listw(analysis_data_nmds_neighbours)    

Moran I statistic standard deviate = 8.0778, p-value = 0.0000000000000003298
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.553407999      -0.003257329       0.004749010 

FDiv HWI

Lat/Long

with_realms_latlong(ggplot(analysis_data, aes(x = latitude, y = longitude, colour = hwi_fdiv_standard)) + geom_point() + standardised_colours_scale + labs(colour = "Standardised response"))

moran.test(analysis_data$hwi_fdiv_standard, nb2listw(analysis_data_neighbours))

    Moran I test under randomisation

data:  analysis_data$hwi_fdiv_standard  
weights: nb2listw(analysis_data_neighbours)    

Moran I statistic standard deviate = 9.7654, p-value < 0.00000000000000022
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.706636768      -0.003257329       0.005284486 

NMDS

with_realms_nmds(ggplot(analysis_data, aes(x = NMDS1, y = NMDS2, colour = hwi_fdiv_standard)) + geom_point() + standardised_colours_scale + labs(colour = "Standardised response"))

moran.test(analysis_data$hwi_fdiv_standard, nb2listw(analysis_data_nmds_neighbours))

    Moran I test under randomisation

data:  analysis_data$hwi_fdiv_standard  
weights: nb2listw(analysis_data_nmds_neighbours)    

Moran I statistic standard deviate = 5.7044, p-value = 0.000000005839
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.389975146      -0.003257329       0.004752081 

FDiv Mass

Lat/Long

with_realms_latlong(ggplot(analysis_data, aes(x = latitude, y = longitude, colour = mass_fdiv_standard)) + geom_point() + standardised_colours_scale + labs(colour = "Standardised response"))

moran.test(analysis_data$mass_fdiv_standard, nb2listw(analysis_data_neighbours))

    Moran I test under randomisation

data:  analysis_data$mass_fdiv_standard  
weights: nb2listw(analysis_data_neighbours)    

Moran I statistic standard deviate = 9.159, p-value < 0.00000000000000022
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.661853721      -0.003257329       0.005273363 

NMDS

with_realms_nmds(ggplot(analysis_data, aes(x = NMDS1, y = NMDS2, colour = mass_fdiv_standard)) + geom_point() + standardised_colours_scale + labs(colour = "Standardised response"))

moran.test(analysis_data$mass_fdiv_standard, nb2listw(analysis_data_nmds_neighbours))

    Moran I test under randomisation

data:  analysis_data$mass_fdiv_standard  
weights: nb2listw(analysis_data_nmds_neighbours)    

Moran I statistic standard deviate = 6.2152, p-value = 0.0000000002563
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.424748274      -0.003257329       0.004742265 

moran.test(analysis_data$mass_fdiv_standard, nb2listw(analysis_data_neighbours))

    Moran I test under randomisation

data:  analysis_data$mass_fdiv_standard  
weights: nb2listw(analysis_data_neighbours)    

Moran I statistic standard deviate = 9.159, p-value < 0.00000000000000022
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.661853721      -0.003257329       0.005273363 

Examine individual metrics

all_explanatories = c(
    'city_avg_ndvi', 'city_avg_elevation', 'city_avg_temp',
    'region_50km_avg_soil_moisture',
    'core_realmAfrotropic', 'core_realmAustralasia', 'core_realmIndomalayan', 'core_realmNearctic', 'core_realmNeotropic', 'core_realmPalearctic',
    'has_introduced_speciesNo introduced species', 'has_introduced_speciesIntroduced species'
)

all_explanatory_names = factor(
   c(
    'Avg. NDVI', 'Avg. Elevation', 'Avg. Temp.',
    'Avg. Soil Moisture',
    'Afrotropic', 'Australasia', 'Indomalayan', 'Nearctic', 'Neotropic', 'Palearctic',
    'Introduced Absent', 'Introduced Present'
  ), ordered = T
)

explanatory_dictionary = data.frame(explanatory = all_explanatories, explanatory_name = all_explanatory_names)
  
with_explanatory_type_labels = function(p) {
  p = p[p$explanatory != '(Intercept)',]
  explanatory_levels = all_explanatories[all_explanatories %in% p$explanatory]
  p$explanatory <- factor(p$explanatory, levels = explanatory_levels)
  
  p$type <- 'Realm'
  p$type[p$explanatory %in% c('city_avg_ndvi', 'city_avg_elevation', 'city_avg_temp')] <- 'City geography'
  p$type[p$explanatory %in% c('region_50km_avg_soil_moisture')] <- 'Regional (50km) geography'
  p$type[p$explanatory %in% c('has_introduced_speciesNo introduced species', 'has_introduced_speciesIntroduced species')] <- 'Introduced species'
  p
}

with_explanatory_names = function(p) {
  p %>% left_join(explanatory_dictionary) %>% arrange(desc(explanatory_name))
}

type_labels = function(p) {
  explanatory_levels = all_explanatories[all_explanatories %in% p$explanatory]
  p$explanatory <- factor(p$explanatory, levels = explanatory_levels)
  
  p$type <- 'Realm'
  p$type[p$explanatory %in% c('city_avg_ndvi', 'city_avg_elevation', 'city_avg_temp', 'city_avg_min_monthly_temp', 'city_avg_max_monthly_temp', 
    'city_avg_monthly_temp', 'city_avg_rainfall', 'city_avg_max_monthly_rainfall', 'city_avg_min_monthly_rainfall', 
    'city_avg_soil_moisture', 'city_max_elev', 'city_min_elev', 'city_elev_range')] <- 'City geography'
  p$type[p$explanatory %in% c('region_50km_avg_ndvi', 'region_50km_avg_elevation', 'region_50km_avg_soil_moisture', 'region_50km_max_elev', 
    'region_50km_min_elev', 'region_50km_elev_range')] <- 'Regional (50km) geography'
   p$type[p$explanatory %in% c('region_20km_avg_ndvi', 'region_20km_avg_elevation', 'region_20km_avg_soil_moisture', 'region_20km_max_elev', 
    'region_20km_min_elev', 'region_20km_elev_range')] <- 'Regional (20km) geography'
  p$type[p$explanatory %in% c('has_introduced_speciesNo introduced species', 'has_introduced_speciesIntroduced species')] <- 'Introduced species'
  p
}

explanatory_labels = c(
  'has_introduced_species'='Introduced species', 
  'has_introduced_speciesNo introduced species'='Introduced absent', 
  'has_introduced_speciesIntroduced species'='Introduced present',
  'city_avg_ndvi'='Average NDVI', 
  'city_avg_elevation'='Average elevation', 
  'city_avg_temp'='Average temperature', 
  'city_avg_min_monthly_temp'='Average minimum monthly temperature', 
  'city_avg_max_monthly_temp'='Average maximum monthly temperature', 
  'city_avg_monthly_temp'='Average monthly temperature', 
  'city_avg_rainfall'='Average rainfall', 
  'city_avg_max_monthly_rainfall'='Average maximum monthly rainfall', 
  'city_avg_min_monthly_rainfall'='Average minimum monthly rainfall', 
  'city_avg_soil_moisture'='Average soil moisture', 
  'city_max_elev'='Maximum elevation', 
  'city_min_elev'='Minimum elevation', 
  'city_elev_range'='Elevation range', 
  'region_20km_avg_ndvi'='Average NDVI', 
  'region_20km_avg_elevation'='Average elevation', 
  'region_20km_avg_soil_moisture'='Average soil moisture', 
  'region_20km_max_elev'='Maximum elevation', 
  'region_20km_min_elev'='Minimum elevation',
  'region_20km_elev_range'='Elevation range',
  'region_50km_avg_ndvi'='Average NDVI',
  'region_50km_avg_elevation'='Average elevation',
  'region_50km_avg_soil_moisture'='Average soil moisture', 
  'region_50km_max_elev'='Maximum elevation',
  'region_50km_min_elev'='Minimum elevation', 
  'region_50km_elev_range'='Elevation range',
  'abs_latitude' = 'Absolute latitude',
  'latitude' = 'Latitude',
  'longitude' = 'Longitude',
  'core_realmAfrotropic' = 'Afrotropical', 
  'core_realmAustralasia' = 'Austaliasian', 
  'core_realmIndomalayan' = 'Indomalayan', 
  'core_realmNearctic' = 'Nearctic', 
  'core_realmNeotropic' = 'Neotropical',
  'core_realmPalearctic' = 'Palearctic',
  'core_realmOceania' = 'Oceanical')

create_formula = function(response_var) {
  as.formula(paste(response_var, '~ core_realm + city_avg_ndvi + city_avg_elevation + city_avg_temp + region_50km_avg_soil_moisture + has_introduced_species'))
}

Helper plot functions

geom_map = function(map_sf, title) {
  norm_mntd_analysis_geo = ggplot() + 
    geom_sf(data = world_map, aes(geometry = geometry)) +
    map_sf +
    standardised_colours_scale +
    labs(colour = 'Standardised\nResponse') +
    theme_bw() +
    theme(legend.position="bottom")
}

Helper Dredge functions

# Taken from MuMIN package
# https://rdrr.io/cran/MuMIn/src/R/averaging.R
# https://rdrr.io/cran/MuMIn/src/R/model.avg.R

.coefarr.avg <-
  function(cfarr, weight, revised.var, full, alpha) {   
    weight <- weight / sum(weight)
    nCoef <- dim(cfarr)[3L]
    if(full) {
      nas <- is.na(cfarr[, 1L, ]) & is.na(cfarr[, 2L, ])
      cfarr[, 1L, ][nas] <- cfarr[, 2L, ][nas] <- 0
      #cfarr[, 1L:2L, ][is.na(cfarr[, 1L:2L, ])] <- 0
      if(!all(is.na(cfarr[, 3L, ])))
        cfarr[ ,3L, ][is.na(cfarr[ , 3L, ])] <- Inf
    }
    
    avgcoef <- array(dim = c(nCoef, 5L),
                     dimnames = list(dimnames(cfarr)[[3L]], c("Estimate",
                                                              "Std. Error", "Adjusted SE", "Lower CI", "Upper CI")))
    for(i in seq_len(nCoef))
      avgcoef[i, ] <- par.avg(cfarr[, 1L, i], cfarr[, 2L, i], weight,
                              df = cfarr[, 3L, i], alpha = alpha, revised.var = revised.var)
    
    avgcoef[is.nan(avgcoef)] <- NA
    return(avgcoef)
  }

.makecoefmat <- function(cf) {
  no.ase <- all(is.na(cf[, 3L]))
  z <- abs(cf[, 1L] / cf[, if(no.ase) 2L else 3L])
  pval <- 2 * pnorm(z, lower.tail = FALSE)
  cbind(cf[, if(no.ase) 1L:2L else 1L:3L, drop = FALSE],
        `z value` = z, `Pr(>|z|)` = zapsmall(pval))
}

# Generate model selections using lmer, dredge, and model.avg
# `forumla` : a two-sided linear formula object describing both the fixed-effects and random-effects part of the model
# `data` : the data frame containing the variables from the formula
# `aic_delta` : the AIC delta to use for selecting models in model average
model_average <- function(formula, data, aic_delta = 20) {
  model <- lm(
    formula,
    data=data
  )
  dredge_result <- dredge(model)
  summary(model.avg(dredge_result, subset = delta < aic_delta))
}

# Create a summary data frame containing the selected variables from a model
# `model_sum` : The model summary output from `model_average`
model_summary <- function(model_sum) {
  .column_name <- function(postfix) {
    postfix
  }
  
  # just return the estimate and p value
  weight <- model_sum$msTable[, 5L]
  
  coefmat.full <- as.data.frame(.makecoefmat(.coefarr.avg(model_sum$coefArray, weight,
                                                          attr(model_sum, "revised.var"), TRUE, 0.05)))
  
  coefmat.subset <-
    as.data.frame(.makecoefmat(.coefarr.avg(model_sum$coefArray, weight,
                                            attr(model_sum, "revised.var"), FALSE, 0.05)))
  
  
  coefmat.subset <- coefmat.subset[-c(1), c(1, 2, 5)]
  names(coefmat.subset) <- c(.column_name("estimate"), .column_name("error"), .column_name("p"))
  coefmat.subset <- tibble::rownames_to_column(coefmat.subset, "explanatory")
  coefmat.subset$model = 'subset'
  
  coefmat.full <- coefmat.full[-c(1), c(1, 2, 5)]
  names(coefmat.full) <- c(.column_name("estimate"), .column_name("error"), .column_name("p"))
  coefmat.full <- tibble::rownames_to_column(coefmat.full, "explanatory")
  coefmat.full$model = 'full'
  
  rbind(coefmat.full, coefmat.subset)
}

plot_dredge_result = function(result_table) {
  p = result_table[result_table$model == 'full',]
  p = type_labels(p)

  ggplot(p, aes(y = explanatory, x = estimate, colour = type)) + 
    geom_line() +
    geom_point() +
    geom_errorbar(aes(xmin=estimate-error, xmax=estimate+error), width=.2,
                   position=position_dodge(0.05)) +
    scale_y_discrete(
      limits = rev(levels(p$explanatory)), 
      labels = explanatory_labels) +
    scale_colour_manual(
      values = c(realm_colour, city_geography_colour, regional_50km_geography_colour, regional_20km_geography_colour, introduced_species_colour), 
      breaks = c('Realm', 'City geography', 'Regional (50km) geography', 'Regional (20km) geography', 'Introduced species')) +
    theme_bw() +
    geom_vline(xintercept=0, linetype="dotted") +
    guides(colour=guide_legend(title="Predictor type")) + xlab('Difference in response from 0\nhabitat filtering (< 0) and competitive interactions (> 0)\n± Standard Error') + ylab('Predictor') +
    theme(legend.justification = "top")
}

GLS Spatial Helpers

gls_method = "ML"

spatial_model = function(formula, correlation) {
  gls(
    formula, 
    data = analysis_data, 
    correlation = correlation, 
    method = gls_method
  )
}

plot_spatial_result = function(model_result) {
  model_summary = summary(model_result)
  result_table = as.data.frame(model_summary$tTable)
  result_table$explanatory = rownames(result_table)
  
  result_table = result_table %>% with_explanatory_type_labels() %>% with_explanatory_names()
  
  ggplot2::ggplot(result_table, ggplot2::aes(y=factor(explanatory_name, level = all_explanatory_names, ordered = T), x=Value, colour = type)) + 
    ggplot2::geom_line() +
    ggplot2::geom_point() +
    ggplot2::geom_errorbar(ggplot2::aes(xmin=Value-Std.Error, xmax=Value+Std.Error), width=.2,
                   position=ggplot2::position_dodge(0.05)) +
    ggplot2::theme_bw() +
    ggplot2::geom_vline(xintercept=0, linetype="dotted") +
    ggplot2::theme(legend.justification = "top") +
    ylab('Predictor') +
    guides(colour=guide_legend(title="Predictor type")) + xlab('Difference in response from 0\nhabitat filtering (< 0) and competitive interactions (> 0)\n± Standard Error') +
    scale_colour_manual(
      values = c(realm_colour, city_geography_colour, regional_50km_geography_colour, introduced_species_colour), 
      breaks = c('Realm', 'City geography', 'Regional (50km) geography', 'Introduced species')) +
    scale_y_discrete(limits = rev(all_explanatory_names[all_explanatory_names %in% result_table$explanatory_name]))
}

Choose best spatial correlation function

AIC(spatial_model(create_formula('mntd_standard'), corLin(form = ~ NMDS1 + NMDS2 + latitude + longitude)))

[1] 662.3448

AIC(spatial_model(create_formula('mntd_standard'), corLin(form = ~ latitude + longitude)))

[1] 668.183

AIC(spatial_model(create_formula('mntd_standard'), corExp(form = ~ NMDS1 + NMDS2 + latitude + longitude)))

[1] 653.2191

AIC(spatial_model(create_formula('mntd_standard'), corExp(form = ~ latitude + longitude)))

[1] 656.8138

AIC(spatial_model(create_formula('mntd_standard'), corGaus(form = ~ NMDS1 + NMDS2 + latitude + longitude)))

[1] 658.2282

AIC(spatial_model(create_formula('mntd_standard'), corGaus(form = ~ latitude + longitude)))

[1] 662.8906

AIC(spatial_model(create_formula('mntd_standard'), corRatio(form = ~ NMDS1 + NMDS2 + latitude + longitude)))

[1] 650.6309

AIC(spatial_model(create_formula('mntd_standard'), corRatio(form = ~ latitude + longitude)))

[1] 655.0259

AIC(spatial_model(create_formula('mntd_standard'), corSpher(form = ~ NMDS1 + NMDS2 + latitude + longitude)))

[1] 660.0803

AIC(spatial_model(create_formula('mntd_standard'), corSpher(form = ~ latitude + longitude)))

[1] 664.3069

MNTD: corRatio with NMDS + lat/long

AIC(spatial_model(create_formula('beak_width_fdiv_standard'), corLin(form = ~ NMDS1 + NMDS2 + latitude + longitude)))

[1] 793.8812

AIC(spatial_model(create_formula('beak_width_fdiv_standard'), corLin(form = ~ latitude + longitude)))

[1] 793.8812

AIC(spatial_model(create_formula('beak_width_fdiv_standard'), corExp(form = ~ NMDS1 + NMDS2 + latitude + longitude)))

[1] 789.9986

AIC(spatial_model(create_formula('beak_width_fdiv_standard'), corExp(form = ~ latitude + longitude)))

[1] 790.7424

AIC(spatial_model(create_formula('beak_width_fdiv_standard'), corGaus(form = ~ NMDS1 + NMDS2 + latitude + longitude)))

[1] 794.1143

AIC(spatial_model(create_formula('beak_width_fdiv_standard'), corGaus(form = ~ latitude + longitude)))

[1] 794.1127

AIC(spatial_model(create_formula('beak_width_fdiv_standard'), corRatio(form = ~ NMDS1 + NMDS2 + latitude + longitude)))

[1] 791.6815

AIC(spatial_model(create_formula('beak_width_fdiv_standard'), corRatio(form = ~ latitude + longitude)))

[1] 791.833

AIC(spatial_model(create_formula('beak_width_fdiv_standard'), corSpher(form = ~ NMDS1 + NMDS2 + latitude + longitude)))

[1] 794.1625

AIC(spatial_model(create_formula('beak_width_fdiv_standard'), corSpher(form = ~ latitude + longitude)))

[1] 794.1625

Beak width: corExp with NMDS + lat/long

AIC(spatial_model(create_formula('hwi_fdiv_standard'), corLin(form = ~ NMDS1 + NMDS2 + latitude + longitude)))

[1] 695.3483

AIC(spatial_model(create_formula('hwi_fdiv_standard'), corLin(form = ~ latitude + longitude)))

[1] 699.9416

AIC(spatial_model(create_formula('hwi_fdiv_standard'), corExp(form = ~ NMDS1 + NMDS2 + latitude + longitude)))

[1] 667.3598

AIC(spatial_model(create_formula('hwi_fdiv_standard'), corExp(form = ~ latitude + longitude)))

[1] 672.0891

AIC(spatial_model(create_formula('hwi_fdiv_standard'), corGaus(form = ~ NMDS1 + NMDS2 + latitude + longitude)))

[1] 705.4322

AIC(spatial_model(create_formula('hwi_fdiv_standard'), corGaus(form = ~ latitude + longitude)))

[1] 709.9593

AIC(spatial_model(create_formula('hwi_fdiv_standard'), corRatio(form = ~ NMDS1 + NMDS2 + latitude + longitude)))

[1] 686.368

AIC(spatial_model(create_formula('hwi_fdiv_standard'), corRatio(form = ~ latitude + longitude)))

[1] 691.6763

AIC(spatial_model(create_formula('hwi_fdiv_standard'), corSpher(form = ~ NMDS1 + NMDS2 + latitude + longitude)))

[1] 690.2192

AIC(spatial_model(create_formula('hwi_fdiv_standard'), corSpher(form = ~ latitude + longitude)))

[1] 696.1892

HWI: corExp with NMDS + lat/long

AIC(spatial_model(create_formula('mass_fdiv_standard'), corLin(form = ~ NMDS1 + NMDS2 + latitude + longitude)))

[1] 777.7153

AIC(spatial_model(create_formula('mass_fdiv_standard'), corLin(form = ~ latitude + longitude)))

[1] 776.4599

AIC(spatial_model(create_formula('mass_fdiv_standard'), corExp(form = ~ NMDS1 + NMDS2 + latitude + longitude)))

[1] 752.6284

AIC(spatial_model(create_formula('mass_fdiv_standard'), corExp(form = ~ latitude + longitude)))

[1] 755.2302

AIC(spatial_model(create_formula('mass_fdiv_standard'), corGaus(form = ~ NMDS1 + NMDS2 + latitude + longitude)))

[1] 786.4347

AIC(spatial_model(create_formula('mass_fdiv_standard'), corGaus(form = ~ latitude + longitude)))

[1] 784.9074

AIC(spatial_model(create_formula('mass_fdiv_standard'), corRatio(form = ~ NMDS1 + NMDS2 + latitude + longitude)))

[1] 766.6933

AIC(spatial_model(create_formula('mass_fdiv_standard'), corRatio(form = ~ latitude + longitude)))

[1] 766.7129

AIC(spatial_model(create_formula('mass_fdiv_standard'), corSpher(form = ~ NMDS1 + NMDS2 + latitude + longitude)))

[1] 780.5551

AIC(spatial_model(create_formula('mass_fdiv_standard'), corSpher(form = ~ latitude + longitude)))

[1] 780.1225

Mass: corExp with NMDS + lat/long

correlation_formula = as.formula('~ NMDS1 + NMDS2 + latitude + longitude')

correlation_function_fdiv = function() {
  corExp(form = correlation_formula)
}

correlation_function_mntd = function() {
  corRatio(form = correlation_formula)
}

MNTD

std_mntd_analysis_geo_plot = geom_map(geom_sf(data = analysis_data, aes(color = mntd_standard, geometry = geometry)), 'MNTD')
std_mntd_analysis_geo_plot

std_mntd_analysis_data = model_data(analysis_data[!is.na(analysis_data$mntd_standard),], 'mntd_standard')
std_mntd_analysis_formula = create_formula('mntd_standard')
std_mntd_analysis_result <- model_average(std_mntd_analysis_formula, std_mntd_analysis_data)

Fixed term is "(Intercept)"

std_mntd_analysis_result_table = model_summary(std_mntd_analysis_result)
std_mntd_analysis_result_table

std_mntd_analysis_pred_plot = plot_dredge_result(std_mntd_analysis_result_table)
std_mntd_analysis_pred_plot

Do the residuals still contain spatial autocorrelation from a fitted lm?

std_mntd_lm = lm(std_mntd_analysis_formula, std_mntd_analysis_data)
moran.test(std_mntd_lm$residuals, nb2listw(analysis_data_neighbours))

    Moran I test under randomisation

data:  std_mntd_lm$residuals  
weights: nb2listw(analysis_data_neighbours)    

Moran I statistic standard deviate = 3.8174, p-value = 0.00006743
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.273851651      -0.003257329       0.005269469 

moran.test(std_mntd_lm$residuals, nb2listw(analysis_data_nmds_neighbours))

    Moran I test under randomisation

data:  std_mntd_lm$residuals  
weights: nb2listw(analysis_data_nmds_neighbours)    

Moran I statistic standard deviate = 2.9313, p-value = 0.001688
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.198532827      -0.003257329       0.004738829 

mntd_spatial_model = spatial_model(std_mntd_analysis_formula, correlation_function_mntd())
moran.test(mntd_spatial_model$residuals, nb2listw(analysis_data_neighbours))

    Moran I test under randomisation

data:  mntd_spatial_model$residuals  
weights: nb2listw(analysis_data_neighbours)    

Moran I statistic standard deviate = 4.0185, p-value = 0.00002929
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.288487485      -0.003257329       0.005270882 

moran.test(mntd_spatial_model$residuals, nb2listw(analysis_data_nmds_neighbours))

    Moran I test under randomisation

data:  mntd_spatial_model$residuals  
weights: nb2listw(analysis_data_nmds_neighbours)    

Moran I statistic standard deviate = 3.1254, p-value = 0.0008877
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.211923502      -0.003257329       0.004740076 

std_mntd_analysis_pred_spatial_plot = plot_spatial_result(mntd_spatial_model)

Joining with `by = join_by(explanatory)`

std_mntd_analysis_pred_spatial_plot

Gape width - FDiv

std_gape_fdiv_analysis_geo_plot = geom_map(geom_sf(data = analysis_data, aes(color = beak_width_fdiv_standard, geometry = geometry)), 'Beak Width FDiv')
std_gape_fdiv_analysis_geo_plot

std_gape_fdiv_analysis_data = model_data(analysis_data[!is.na(analysis_data$beak_width_fdiv_standard),], 'beak_width_fdiv_standard')
std_gape_fdiv_analysis_formula = create_formula('beak_width_fdiv_standard')
std_gape_fdiv_analysis_result <- model_average(std_gape_fdiv_analysis_formula, std_gape_fdiv_analysis_data)

Fixed term is "(Intercept)"

std_gape_fdiv_analysis_result_table = model_summary(std_gape_fdiv_analysis_result)
std_gape_fdiv_analysis_result_table

std_gape_fdiv_analysis_pred_plot = plot_dredge_result(std_gape_fdiv_analysis_result_table)
std_gape_fdiv_analysis_pred_plot

Do the residuals still contain spatial autocorrelation from a fitted lm?

std_gape_fdiv_lm = lm(std_gape_fdiv_analysis_formula, std_gape_fdiv_analysis_data)
moran.test(std_gape_fdiv_lm$residuals, nb2listw(analysis_data_neighbours))

    Moran I test under randomisation

data:  std_gape_fdiv_lm$residuals  
weights: nb2listw(analysis_data_neighbours)    

Moran I statistic standard deviate = 3.55, p-value = 0.0001926
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.254212474      -0.003257329       0.005260213 

moran.test(std_gape_fdiv_lm$residuals, nb2listw(analysis_data_nmds_neighbours))

    Moran I test under randomisation

data:  std_gape_fdiv_lm$residuals  
weights: nb2listw(analysis_data_nmds_neighbours)    

Moran I statistic standard deviate = 1.6631, p-value = 0.04815
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.111128951      -0.003257329       0.004730661 

std_gape_fdiv_spatial_model = spatial_model(std_gape_fdiv_analysis_formula, correlation_function_fdiv())
moran.test(std_gape_fdiv_spatial_model$residuals, nb2listw(analysis_data_neighbours))

    Moran I test under randomisation

data:  std_gape_fdiv_spatial_model$residuals  
weights: nb2listw(analysis_data_neighbours)    

Moran I statistic standard deviate = 3.6108, p-value = 0.0001526
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.258639556      -0.003257329       0.005260832 

moran.test(std_gape_fdiv_spatial_model$residuals, nb2listw(analysis_data_nmds_neighbours))

    Moran I test under randomisation

data:  std_gape_fdiv_spatial_model$residuals  
weights: nb2listw(analysis_data_nmds_neighbours)    

Moran I statistic standard deviate = 1.8619, p-value = 0.03131
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.124813501      -0.003257329       0.004731207 

std_gape_fdiv_analysis_pred_spatial_plot = plot_spatial_result(std_gape_fdiv_spatial_model)

Joining with `by = join_by(explanatory)`

std_gape_fdiv_analysis_pred_spatial_plot

HWI - FDiv

std_hwi_fdiv_analysis_geo_plot = geom_map(geom_sf(data = analysis_data, aes(color = hwi_fdiv_standard, geometry = geometry)), 'HWI FDiv')
std_hwi_fdiv_analysis_geo_plot

std_hwi_fdiv_analysis_data = model_data(analysis_data[!is.na(analysis_data$hwi_fdiv_standard),], 'hwi_fdiv_standard')
std_hwi_fdiv_analysis_formula = create_formula('hwi_fdiv_standard')
std_hwi_fdiv_analysis_result <- model_average(std_hwi_fdiv_analysis_formula, std_hwi_fdiv_analysis_data)

Fixed term is "(Intercept)"

std_hwi_fdiv_analysis_result_table = model_summary(std_hwi_fdiv_analysis_result)
std_hwi_fdiv_analysis_result_table

std_hwi_fdiv_analysis_pred_plot = plot_dredge_result(std_hwi_fdiv_analysis_result_table)
std_hwi_fdiv_analysis_pred_plot

std_hwi_fdiv_lm = lm(std_hwi_fdiv_analysis_formula, std_hwi_fdiv_analysis_data)
moran.test(std_hwi_fdiv_lm$residuals, nb2listw(analysis_data_neighbours))

    Moran I test under randomisation

data:  std_hwi_fdiv_lm$residuals  
weights: nb2listw(analysis_data_neighbours)    

Moran I statistic standard deviate = 6.197, p-value = 0.0000000002877
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.445978540      -0.003257329       0.005255093 

moran.test(std_hwi_fdiv_lm$residuals, nb2listw(analysis_data_nmds_neighbours))

    Moran I test under randomisation

data:  std_hwi_fdiv_lm$residuals  
weights: nb2listw(analysis_data_nmds_neighbours)    

Moran I statistic standard deviate = 3.7792, p-value = 0.00007868
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.256549224      -0.003257329       0.004726143 

std_hwi_fdiv_spatial_model = spatial_model(std_hwi_fdiv_analysis_formula, correlation_function_fdiv())
moran.test(std_hwi_fdiv_spatial_model$residuals, nb2listw(analysis_data_neighbours))

    Moran I test under randomisation

data:  std_hwi_fdiv_spatial_model$residuals  
weights: nb2listw(analysis_data_neighbours)    

Moran I statistic standard deviate = 7.2343, p-value = 0.000000000000234
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.521468855      -0.003257329       0.005261042 

moran.test(std_hwi_fdiv_spatial_model$residuals, nb2listw(analysis_data_nmds_neighbours))

    Moran I test under randomisation

data:  std_hwi_fdiv_spatial_model$residuals  
weights: nb2listw(analysis_data_nmds_neighbours)    

Moran I statistic standard deviate = 3.3574, p-value = 0.0003935
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.227679202      -0.003257329       0.004731393 

std_hwi_fdiv_analysis_pred_spatial_plot = plot_spatial_result(std_hwi_fdiv_spatial_model)

Joining with `by = join_by(explanatory)`

std_hwi_fdiv_analysis_pred_spatial_plot

Mass - FDiv

std_mass_fdiv_analysis_geo_plot = geom_map(geom_sf(data = analysis_data, aes(color = mass_fdiv_standard, geometry = geometry)), 'Mass FDiv')
std_mass_fdiv_analysis_geo_plot

std_mass_fdiv_analysis_data = model_data(analysis_data[!is.na(analysis_data$mass_fdiv_standard),], 'mass_fdiv_standard')
std_mass_fdiv_analysis_formula = create_formula('mass_fdiv_standard')
std_mass_fdiv_analysis_result <- model_average(std_mass_fdiv_analysis_formula, std_mass_fdiv_analysis_data)

Fixed term is "(Intercept)"

std_mass_fdiv_analysis_result_table = model_summary(std_mass_fdiv_analysis_result)
std_mass_fdiv_analysis_result_table

std_mass_fdiv_analysis_pred_plot = plot_dredge_result(std_mass_fdiv_analysis_result_table)
std_mass_fdiv_analysis_pred_plot

std_mass_fdiv_lm = lm(std_mass_fdiv_analysis_formula, std_mass_fdiv_analysis_data)
moran.test(std_mass_fdiv_lm$residuals, nb2listw(analysis_data_neighbours))

    Moran I test under randomisation

data:  std_mass_fdiv_lm$residuals  
weights: nb2listw(analysis_data_neighbours)    

Moran I statistic standard deviate = 7.4464, p-value = 0.00000000000004797
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.536458696      -0.003257329       0.005253399 

moran.test(std_mass_fdiv_lm$residuals, nb2listw(analysis_data_nmds_neighbours))

    Moran I test under randomisation

data:  std_mass_fdiv_lm$residuals  
weights: nb2listw(analysis_data_nmds_neighbours)    

Moran I statistic standard deviate = 4.7703, p-value = 0.0000009196
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.324636672      -0.003257329       0.004724648 

std_mass_fdiv_spatial_model = spatial_model(std_mass_fdiv_analysis_formula, correlation_function_fdiv())
moran.test(std_mass_fdiv_spatial_model$residuals, nb2listw(analysis_data_neighbours))

    Moran I test under randomisation

data:  std_mass_fdiv_spatial_model$residuals  
weights: nb2listw(analysis_data_neighbours)    

Moran I statistic standard deviate = 7.9834, p-value = 0.0000000000000007117
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.575429652      -0.003257329       0.005254232 

moran.test(std_mass_fdiv_spatial_model$residuals, nb2listw(analysis_data_nmds_neighbours))

    Moran I test under randomisation

data:  std_mass_fdiv_spatial_model$residuals  
weights: nb2listw(analysis_data_nmds_neighbours)    

Moran I statistic standard deviate = 5.2794, p-value = 0.00000006481
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.359654511      -0.003257329       0.004725382 

std_mass_fdiv_analysis_pred_spatial_plot = plot_spatial_result(std_mass_fdiv_spatial_model)

Joining with `by = join_by(explanatory)`

std_mass_fdiv_analysis_pred_spatial_plot

Create plot of differences in process response

pred_legend <- ggpubr::get_legend(
  # create some space to the left of the legend
  std_hwi_fdiv_analysis_pred_plot + theme(legend.box.margin = margin(0, 0, 0, 0)) + guides(colour=guide_legend(ncol=2)) + labs(color = "Predictor type")
)

`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?

geo_legend <- ggpubr::get_legend(
  # create some space to the left of the legend
  std_mass_fdiv_analysis_geo_plot + theme(legend.box.margin = margin(-80, 0, 0, 12), legend.title.position = "top", legend.key.width = unit(10, 'mm')) + labs(color = "Standardised response")
)

legend = plot_grid(
  geo_legend,
  pred_legend, 
  nrow = 1
)
legend

plot_grid(
  plot_grid(
    std_mntd_analysis_geo_plot + theme(legend.position="none"), 
    std_mntd_analysis_pred_plot + theme(legend.position="none") + scale_x_continuous(name = '', limits = c(-3, 3)) + ylab(''), 
    nrow = 1
  ) + draw_label("MNTD", size = 16, angle = 90, x = 0.01, y = 0.5),
  plot_grid(
    std_gape_fdiv_analysis_geo_plot + theme(legend.position="none"), 
    std_gape_fdiv_analysis_pred_plot + theme(legend.position="none") + scale_x_continuous(name = '', limits = c(-3, 3)) + ylab(''), 
    nrow = 1
  ) + draw_label("Beak Width", size = 16, angle = 90, x = 0.01, y = 0.5),
  plot_grid(
    std_hwi_fdiv_analysis_geo_plot + theme(legend.position="none"), 
    std_hwi_fdiv_analysis_pred_plot + theme(legend.position="none") + scale_x_continuous(name = '', limits = c(-3, 3)) + ylab(''), 
    nrow = 1
  ) + draw_label("HWI", size = 16, angle = 90, x = 0.01, y = 0.5),
  plot_grid(
    std_mass_fdiv_analysis_geo_plot + theme(legend.position="none"), 
    std_mass_fdiv_analysis_pred_plot + theme(legend.position="none") + scale_x_continuous(name = '', limits = c(-3, 3)) + ylab(''), 
    nrow = 1
  ) + draw_label("Mass", size = 16, angle = 90, x = 0.01, y = 0.5), 
  legend,
  nrow = 5
)

`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?

ggsave(filename(FIGURES_OUTPUT_DIR, 'process_response.jpg'), width = 3000, height = 3200, units = 'px')

Compare to spatial model

plot_grid(
  plot_grid(
    ggdraw() + 
    draw_label(
      "Spatial Model",
      fontface = 'bold',
      x = 0,
      hjust = 0
    ),
    ggdraw() + 
    draw_label(
      "Dredge Result",
      fontface = 'bold',
      x = 0,
      hjust = 0
    )
  ),
  plot_grid(
    std_mntd_analysis_pred_spatial_plot + theme(legend.position="none") + scale_x_continuous(name = '', limits = c(-3, 3)) + ylab(''), 
    std_mntd_analysis_pred_plot + theme(legend.position="none") + scale_x_continuous(name = '', limits = c(-3, 3)) + ylab(''), 
    nrow = 1
  ) + draw_label("MNTD", size = 16, angle = 90, x = 0.01, y = 0.5),
  plot_grid(
    std_gape_fdiv_analysis_pred_spatial_plot + theme(legend.position="none") + scale_x_continuous(name = '', limits = c(-3, 3)) + ylab(''),  
    std_gape_fdiv_analysis_pred_plot + theme(legend.position="none") + scale_x_continuous(name = '', limits = c(-3, 3)) + ylab(''), 
    nrow = 1
  ) + draw_label("Beak Width", size = 16, angle = 90, x = 0.01, y = 0.5),
  plot_grid(
    std_hwi_fdiv_analysis_pred_spatial_plot + theme(legend.position="none") + scale_x_continuous(name = '', limits = c(-3, 3)) + ylab(''), 
    std_hwi_fdiv_analysis_pred_plot + theme(legend.position="none") + scale_x_continuous(name = '', limits = c(-3, 3)) + ylab(''), 
    nrow = 1
  ) + draw_label("HWI", size = 16, angle = 90, x = 0.01, y = 0.5),
  plot_grid(
    std_mass_fdiv_analysis_pred_spatial_plot + theme(legend.position="none") + scale_x_continuous(name = '', limits = c(-3, 3)) + ylab(''), 
    std_mass_fdiv_analysis_pred_plot + theme(legend.position="none") + scale_x_continuous(name = '', limits = c(-3, 3)) + ylab(''), 
    nrow = 1
  ) + draw_label("Mass", size = 16, angle = 90, x = 0.01, y = 0.5), 
  nrow = 5
)

`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?

ggsave(filename(FIGURES_OUTPUT_DIR, 'process_response_vs_spatial.jpg'), width = 3000, height = 3200, units = 'px')

Compare metrics against each other

ggplot(analysis_data, aes(x = beak_width_fdiv_standard, y = mntd_standard, colour = core_realm)) + 
  geom_point() +
  ylab("MNTD") + 
  xlab("Beak Width FDiv") +
  theme_bw() + labs(color = "Realm")

ggplot(analysis_data, aes(x = hwi_fdiv_standard, y = mntd_standard, colour = core_realm)) + 
  geom_point() +
  ylab("MNTD") + 
  xlab("HWI FDiv") +
  theme_bw() + labs(color = "Realm")

ggplot(analysis_data, aes(x = hwi_fdiv_standard, y = beak_width_fdiv_standard, colour = core_realm)) + 
  geom_point() +
  ylab("Beak Width FDiv") + 
  xlab("HWI FDiv") +
  theme_bw() + labs(color = "Realm")

mntd_fdiv_analysis = analysis_data %>% 
  dplyr::select(city_id,  mntd_standard, hwi_fdiv_standard, beak_width_fdiv_standard, mass_fdiv_standard) %>%
  left_join(community_summary) %>%
  mutate(urban_pool_perc = urban_pool_size * 100 / regional_pool_size)

Joining with `by = join_by(city_id)`

mntd_fdiv_analysis

ggpairs(mntd_fdiv_analysis %>% dplyr::select(mntd_standard, hwi_fdiv_standard, beak_width_fdiv_standard, mass_fdiv_standard, regional_pool_size, urban_pool_size, urban_pool_perc), columnLabels = c('MNTD', 'HWI FD', 'Bk FD', 'Mss FD', 'Region Rich.', 'Urban Rich.', '% Urban'))
ggsave(filename(FIGURES_OUTPUT_DIR, 'appendix_standarised_correlation.jpg'))

Saving 7.29 x 4.51 in image